//
// File: calc_TimeFreq_Ar.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 21-Mar-2025 17:03:22
//

// Include Files
#include "calc_TimeFreq_Ar.h"
#include "calcAR.h"
#include "calcAreas.h"
#include "calc_TimeFreq_Ar_data.h"
#include "calc_TimeFreq_Ar_initialize.h"
#include "calc_TimeFreq_Ar_types.h"
#include "colon.h"
#include "hamming.h"
#include "interp1.h"
#include "mean.h"
#include "pburg.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include <cmath>

// Function Definitions
//
// Arguments    : const coder::array<double, 1U> &t
//                const coder::array<double, 1U> &y
//                double fs
//                double nfft
//                double AR_order
//                double winSize
//                double overlap
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct0_T *ar_out
// Return Type  : void
//
void calc_TimeFreq_Ar(const coder::array<double, 1U> &t,
                      const coder::array<double, 1U> &y, double fs, double nfft,
                      double AR_order, double winSize, double overlap,
                      const double VLF[2], const double LF[2],
                      const double HF[2], struct0_T *ar_out)
{
  coder::array<double, 2U> iStart;
  coder::array<double, 2U> idx;
  coder::array<double, 2U> tint;
  coder::array<double, 1U> b_y;
  coder::array<double, 1U> f;
  coder::array<double, 1U> y2;
  coder::array<boolean_T, 1U> x;
  struct2_T expl_temp;
  double b_winSize;
  double bsum;
  double xtmp;
  int i;
  int i1;
  int k;
  int loop_ub;
  int nblocks;
  int nd2;
  int nx;
  boolean_T b;
  if (!isInitialized_calc_TimeFreq_Ar) {
    calc_TimeFreq_Ar_initialize();
  }
  // calAR - Calculates PSD using windowed Burg method.
  //
  // Inputs:
  // Outputs:
  b_winSize = winSize * fs;
  //  (samples)
  // resample
  xtmp = 1.0 / fs;
  if (std::isnan(t[0]) || std::isnan(xtmp) || std::isnan(t[t.size(0) - 1])) {
    tint.set_size(1, 1);
    tint[0] = rtNaN;
  } else if ((xtmp == 0.0) || ((t[0] < t[t.size(0) - 1]) && (xtmp < 0.0)) ||
             ((t[t.size(0) - 1] < t[0]) && (xtmp > 0.0))) {
    tint.set_size(1, 0);
  } else if ((std::isinf(t[0]) || std::isinf(t[t.size(0) - 1])) &&
             (std::isinf(xtmp) || (t[0] == t[t.size(0) - 1]))) {
    tint.set_size(1, 1);
    tint[0] = rtNaN;
  } else if (std::isinf(xtmp)) {
    tint.set_size(1, 1);
    tint[0] = t[0];
  } else if ((std::floor(t[0]) == t[0]) && (std::floor(xtmp) == xtmp)) {
    tint.set_size(1, static_cast<int>((t[t.size(0) - 1] - t[0]) / xtmp) + 1);
    loop_ub = static_cast<int>((t[t.size(0) - 1] - t[0]) / xtmp);
    for (i = 0; i <= loop_ub; i++) {
      tint[i] = t[0] + xtmp * static_cast<double>(i);
    }
  } else {
    coder::eml_float_colon(t[0], xtmp, t[t.size(0) - 1], tint);
  }
  // time values for interp.
  b_y.set_size(y.size(0));
  loop_ub = y.size(0);
  for (i = 0; i < loop_ub; i++) {
    b_y[i] = y[i];
  }
  ar_out->global_PSD.set_size(t.size(0));
  loop_ub = t.size(0);
  for (i = 0; i < loop_ub; i++) {
    ar_out->global_PSD[i] = t[i];
  }
  nx = t.size(0) - 1;
  iStart.set_size(1, tint.size(1));
  loop_ub = tint.size(1);
  for (i = 0; i < loop_ub; i++) {
    iStart[i] = 0.0;
  }
  b = (tint.size(1) == 0);
  if (!b) {
    k = 0;
    int exitg1;
    do {
      exitg1 = 0;
      if (k <= nx) {
        if (std::isnan(t[k])) {
          exitg1 = 1;
        } else {
          k++;
        }
      } else {
        if (t[1] < t[0]) {
          i = (nx + 1) >> 1;
          for (nd2 = 0; nd2 < i; nd2++) {
            xtmp = ar_out->global_PSD[nd2];
            i1 = nx - nd2;
            ar_out->global_PSD[nd2] = ar_out->global_PSD[i1];
            ar_out->global_PSD[i1] = xtmp;
          }
          if (y.size(0) > 1) {
            nx = y.size(0) - 1;
            nd2 = y.size(0) >> 1;
            for (k = 0; k < nd2; k++) {
              xtmp = b_y[k];
              nblocks = nx - k;
              b_y[k] = b_y[nblocks];
              b_y[nblocks] = xtmp;
            }
          }
        }
        coder::interp1SplineOrPCHIP(b_y, tint, iStart, ar_out->global_PSD);
        exitg1 = 1;
      }
    } while (exitg1 == 0);
  }
  b_y.set_size(iStart.size(1));
  loop_ub = iStart.size(1);
  for (i = 0; i < loop_ub; i++) {
    b_y[i] = iStart[i];
  }
  // cubic spline interpolation
  // get limits of windows
  if (tint[tint.size(1) - 1] >= b_winSize) {
    // windowArray - Seperates the input array s into smaller segments/windows
    // with lengths of winSize. winSize can be either in units of time or
    // samples. If winSize is given in units of time the first column of the
    // input array must containg time and method should be set to a value of 1.
    //
    // INPUTS:
    //    s: input array OR length of array. If outputing segs or using the time
    //        methods, the 1st dim must contain time values.
    //    winSize: size of segments to return
    //    method: method of segmenting array (0=by samples, 1=by time)
    //    overlap: amount that segments/windows overlap.
    //    flagCenters: decides whether to output center of windows instead of
    //    beginning
    //
    // OUTPUTS:
    //    idx: 2dim array that contains the start and end of each segment/window
    //    centers: array containg locations of center of windows (useful when
    //    using the time method)
    //
    // EXAMPLES:
    //    1. Seperate input using window size of 100 seconds and overlap
    //    of 50 seconds.
    //    idx=slidingWindow(input,100,50,1);
    //
    //    2. Seperate input using window size of 128 samples and overlap
    //    of 64 samples.
    //    idx=slidingWindow(input,128,64,0);
    //
    //    3. Seperate input using window size of 128 samples and overlap
    //    of 64 samples. Also outputs a cell array containing the actual
    //    segments [idx, segs]=slidingWindow(input,128,64,0,1);
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
    //
    //  This file is part of HRVAS
    //
    //  HRVAS is free software: you can redistribute it and/or modify
    //  it under the terms of the GNU General Public License as published by
    //  the Free Software Foundation, either version 3 of the License, or
    //  (at your option) any later version.
    //
    //  HRVAS is distributed in the hope that it will be useful,
    //  but WITHOUT ANY WARRANTY; without even the implied warranty of
    //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    //  GNU General Public License for more details.
    //
    //  You should have received a copy of the GNU General Public License
    //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // Check inputs
    // includeSeg=0;
    // check inputs
    //  if (m==1) && (n==1) && (method==1)
    //      error('slidingWindow: Invalid input array for this windowing
    //      method.') return;
    //  end
    //
    //  if (m==1) && (n==1)
    //      N=s;
    //  %     if includeSeg
    //  %         warning('SlidingWindow: Cannot output segs if input is
    //  scalar'); %     end
    //      includeSeg=false;
    //  else
    //      N=length(s);
    //  end
    // segment/window data based on samples
    xtmp = b_winSize - overlap * fs;
    bsum = (static_cast<double>(tint.size(1)) - b_winSize) + 1.0;
    if (std::isnan(xtmp)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if ((xtmp == 0.0) || ((bsum > 1.0) && (xtmp < 0.0)) ||
               ((bsum < 1.0) && (xtmp > 0.0))) {
      iStart.set_size(1, 0);
    } else if (std::isinf(bsum) && std::isinf(xtmp)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if (std::isinf(xtmp)) {
      iStart.set_size(1, 1);
      iStart[0] = 1.0;
    } else if (std::floor(xtmp) == xtmp) {
      loop_ub = static_cast<int>((bsum - 1.0) / xtmp);
      iStart.set_size(1, loop_ub + 1);
      for (i = 0; i <= loop_ub; i++) {
        iStart[i] = xtmp * static_cast<double>(i) + 1.0;
      }
    } else {
      coder::eml_float_colon(xtmp, bsum, iStart);
    }
    idx.set_size(iStart.size(1), 2);
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i] = iStart[i];
    }
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i + idx.size(0)] = iStart[i] + (b_winSize - 1.0);
    }
    //  if includeSeg
    //      segs=cell(size(idx,1));
    //      for i=1:size(idx,1)
    //          segs{i}=s(idx(i,1):idx(i,2));%code here for outputing
    //      end
    //  end
    // (sample #)
  } else {
    idx.set_size(1, 2);
    idx[0] = 1.0;
    idx[1] = tint.size(1);
  }
  xtmp = std::round(b_winSize / 2.0);
  loop_ub = idx.size(0);
  ar_out->T.set_size(1, idx.size(0));
  for (i = 0; i < loop_ub; i++) {
    ar_out->T[i] = tint[static_cast<int>(idx[i] + xtmp) - 1];
  }
  // calculate center time of window (s)
  // used for plotting
  // preallocate memory
  // number of PSD/windows
  ar_out->psd.set_size(static_cast<int>(nfft), idx.size(0));
  ar_out->f.set_size(static_cast<int>(nfft), idx.size(0));
  loop_ub = static_cast<int>(nfft) * idx.size(0);
  for (i = 0; i < loop_ub; i++) {
    ar_out->f[i] = 0.0;
  }
  // Calculate PSD
  i = idx.size(0);
  for (int b_i{0}; b_i < i; b_i++) {
    // Prepare y2 and t2
    xtmp = idx[b_i];
    bsum = idx[b_i + idx.size(0)];
    if (xtmp > bsum) {
      i1 = 0;
      nd2 = 0;
    } else {
      i1 = static_cast<int>(xtmp) - 1;
      nd2 = static_cast<int>(bsum);
    }
    //          t2=tint(idx(i,1):idx(i,2));
    // remove linear trend
    //         y=detrend(y,'linear');
    loop_ub = nd2 - i1;
    if (loop_ub == 0) {
      xtmp = 0.0;
    } else {
      int lastBlockLength;
      if (loop_ub <= 1024) {
        nd2 = loop_ub;
        lastBlockLength = 0;
        nblocks = 1;
      } else {
        nd2 = 1024;
        nblocks = loop_ub / 1024;
        lastBlockLength = loop_ub - (nblocks << 10);
        if (lastBlockLength > 0) {
          nblocks++;
        } else {
          lastBlockLength = 1024;
        }
      }
      xtmp = b_y[i1];
      for (k = 2; k <= nd2; k++) {
        xtmp += b_y[(i1 + k) - 1];
      }
      for (int ib{2}; ib <= nblocks; ib++) {
        nd2 = i1 + ((ib - 1) << 10);
        bsum = b_y[nd2];
        if (ib == nblocks) {
          nx = lastBlockLength;
        } else {
          nx = 1024;
        }
        for (k = 2; k <= nx; k++) {
          bsum += b_y[(nd2 + k) - 1];
        }
        xtmp += bsum;
      }
    }
    xtmp /= static_cast<double>(loop_ub - 1) + 1.0;
    y2.set_size(loop_ub);
    for (nd2 = 0; nd2 < loop_ub; nd2++) {
      y2[nd2] = b_y[i1 + nd2] - xtmp;
    }
    // remove mean
    coder::hamming(static_cast<double>(y2.size(0)), ar_out->global_PSD);
    if (y2.size(0) == ar_out->global_PSD.size(0)) {
      loop_ub = y2.size(0);
      for (i1 = 0; i1 < loop_ub; i1++) {
        y2[i1] = y2[i1] * ar_out->global_PSD[i1];
      }
    } else {
      binary_expand_op(y2, ar_out);
    }
    // hamming window
    // Calculate PSD
    coder::pburg(y2, AR_order, nfft * 2.0 - 1.0, fs, ar_out->global_PSD, f);
    loop_ub = f.size(0);
    ar_out->f.set_size(f.size(0), 1);
    for (i1 = 0; i1 < loop_ub; i1++) {
      ar_out->f[i1] = f[i1];
    }
    loop_ub = ar_out->global_PSD.size(0);
    for (i1 = 0; i1 < loop_ub; i1++) {
      ar_out->psd[i1 + ar_out->psd.size(0) * b_i] = ar_out->global_PSD[i1];
    }
  }
  //  calcAreas - Calulates areas/energy under the PSD curve within the freq
  //  bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
  //  percentage, and normalized units. Also returns LF/HF ratio.
  //
  //  Inputs:
  //    PSD: PSD vector
  //    F: Freq vector
  //    VLF, LF, HF: array containing VLF, LF, and HF freq limits
  //    flagVLF: flag to decide whether to calculate VLF hrv
  //  Output:
  //
  //  Usage:
  //
  //
  //  Ref: Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
  ar_out->hrv.aVLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.aVLF[i] = 0.0;
  }
  ar_out->hrv.aLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.aLF[i] = 0.0;
  }
  ar_out->hrv.aHF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.aHF[i] = 0.0;
  }
  ar_out->hrv.aTotal.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.aTotal[i] = 0.0;
  }
  ar_out->hrv.pVLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.pVLF[i] = 0.0;
  }
  ar_out->hrv.pLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.pLF[i] = 0.0;
  }
  ar_out->hrv.pHF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.pHF[i] = 0.0;
  }
  ar_out->hrv.nLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.nLF[i] = 0.0;
  }
  ar_out->hrv.nHF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.nHF[i] = 0.0;
  }
  ar_out->hrv.LFHF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.LFHF[i] = 0.0;
  }
  ar_out->hrv.peakVLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.peakVLF[i] = 0.0;
  }
  ar_out->hrv.peakLF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.peakLF[i] = 0.0;
  }
  ar_out->hrv.peakHF.set_size(ar_out->psd.size(1));
  loop_ub = ar_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->hrv.peakHF[i] = 0.0;
  }
  i = ar_out->psd.size(1);
  for (nd2 = 0; nd2 < i; nd2++) {
    loop_ub = ar_out->psd.size(0);
    b_y.set_size(ar_out->psd.size(0));
    for (i1 = 0; i1 < loop_ub; i1++) {
      b_y[i1] = ar_out->psd[i1 + ar_out->psd.size(0) * nd2];
    }
    calcAreas(ar_out->f, b_y, VLF, LF, HF, &expl_temp);
    // create output structure
    ar_out->hrv.aVLF[nd2] = expl_temp.aVLF;
    ar_out->hrv.aLF[nd2] = expl_temp.aLF;
    ar_out->hrv.aHF[nd2] = expl_temp.aHF;
    ar_out->hrv.aTotal[nd2] = expl_temp.aTotal;
    ar_out->hrv.pVLF[nd2] = expl_temp.pVLF;
    ar_out->hrv.pLF[nd2] = expl_temp.pLF;
    ar_out->hrv.pHF[nd2] = expl_temp.pHF;
    ar_out->hrv.nLF[nd2] = expl_temp.nLF;
    ar_out->hrv.nHF[nd2] = expl_temp.nHF;
    ar_out->hrv.LFHF[nd2] = expl_temp.LFHF;
    ar_out->hrv.peakVLF[nd2] = expl_temp.peakVLF;
    ar_out->hrv.peakLF[nd2] = expl_temp.peakLF;
    ar_out->hrv.peakHF[nd2] = expl_temp.peakHF;
  }
  x.set_size(ar_out->hrv.LFHF.size(0));
  loop_ub = ar_out->hrv.LFHF.size(0);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (ar_out->hrv.LFHF[i] > 1.0);
  }
  nx = x.size(0);
  if (x.size(0) == 0) {
    nblocks = 0;
  } else {
    nblocks = x[0];
    for (k = 2; k <= nx; k++) {
      nblocks += x[k - 1];
    }
  }
  x.set_size(ar_out->hrv.LFHF.size(0));
  loop_ub = ar_out->hrv.LFHF.size(0);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (ar_out->hrv.LFHF[i] <= 1.0);
  }
  nx = x.size(0);
  if (x.size(0) == 0) {
    nd2 = 0;
  } else {
    nd2 = x[0];
    for (k = 2; k <= nx; k++) {
      nd2 += x[k - 1];
    }
  }
  // global psd
  coder::mean(ar_out->psd, ar_out->global_PSD);
  calcAreas(ar_out->f, ar_out->global_PSD, VLF, LF, HF, &ar_out->global_hrv);
  ar_out->hrv.rLFHF = static_cast<double>(nblocks) / static_cast<double>(nd2);
  ar_out->global_f.set_size(ar_out->f.size(0), ar_out->f.size(1));
  loop_ub = ar_out->f.size(0) * ar_out->f.size(1);
  for (i = 0; i < loop_ub; i++) {
    ar_out->global_f[i] = ar_out->f[i];
  }
}

//
// File trailer for calc_TimeFreq_Ar.cpp
//
// [EOF]
//
